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Abstract 

The hydrodynamics of viscoelastic materials (for example polymer melts and solu- 
tions) presents interesting and complex phenomena, for example instabilities and 
turbulent flow at very low Reynolds numbers due to normal stress effects and the ex- 
istence of a finite stress relaxation time. This present work is motivated by renewed 
interest in instabilities in polymer flow. The majority of currently used numerical 
methods discretize a constitutive equation on a grid with finite difference or similar 
methods. We present work in progress in which we simulate viscoelastic flow with 
dissipative particle dynamics. The advantage of this approach is that many of the 
numerical instabilities of conventional methods can be avoided, and that the model 
gives clear physical insight into the origins of many viscoelastic flow instabilities. 
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1 Introduction 



One of the most important properties of viscoelastic fluids is that shear flow 
affects not only the off-diagonal (shear) component of the stress tensor, but 
also the diagonal elements: they change with respect to each other [1,2]. In 
the plane Couette flow geometry of Fig. 1, where the viscosity rj is 

<Txz = vi (i) 
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Fig. 1. The plane Couette flow geometry of our simulations. 



(cr a i3 is the stress tensor and 7 = dv x /dz is the shear rate), the first normal 
stress difference N\ is defined as 

<Txx-<Tzz = ^1(7) = ^17 2 , (2) 

where the first normal stress coefficient ^1 is finite for small shear rate 7. 

Many of the unusual flow phenomena of viscoelastic fluids can be traced back 
to this nonzero normal stress difference. The extra forces generated by the vis- 
coelastic stresses often destabilize the flow, leading to instabilities. This work 
has been motivated by the realization that many key issues regarding vis- 
coelastic instabilities are unresolved. Eg., it was suggested recently that while 
the viscoelastic Poiseuille flow might be linearly stable, it could be nonlinearly 
unstable, and this might be a route to melt-fracture type behavior [3,4]. These 
viscoelastic instabilities occur even at vanishing Reynolds numbers, and are 
driven not by kinetic forces but by elastic forces. The control parameter is 
the Weissenberg number Wi, the ratio of the characteristic relaxation time of 
the fluid (r) and the characteristic time of the flow: for a shear flow it is 77" . 
As Wi is increased, the instabilities lead to a complex non- stationary flow, 
often called viscoelastic turbulence [5,6], based on similarities to the Reynolds 
number driven turbulence. 

Clearly numerical methods are in great need to understand these complex 
phenomena. However, the classical rheological engineering approach, based 
on finite volume or finite difference discretization of the Navier-Stokes equa- 
tion together with the viscoelastic constitutive equations, runs into numerical 
instabilities in what is commonly referred to as "high Weissenberg number 
problem" [7]. Our goal is to overcome this barrier by turning to alternative, 
discrete methods of fluid dynamics. 

In this paper we describe an extension of one of the successful discrete meth- 
ods, the dissipative particle dynamics, to viscoelastic fluids. In this approach 
we don't start from a closed set of equations, instead use kinetic considera- 
tions. This is work in progress, and we show the first validation results. 



2 Dissipative particle dynamics 



In the method of dissipative particle dynamics (DPD) [8,9] the fluid is rep- 
resented with particles, each one corresponding to a macroscopic blob of the 
fluid. The particles interact with each other via finite range pairwise forces. 
The force exerted on particle % of a Newtonian fluid can be written as a sum 
of conservative, dissipative and random contributions: 



f, = E ( f r s + *S iss + f r d ) 



(3) 



The conservative part of the force is a soft repulsion: 



.peons 

ij ~ S 



o, 



(4) 



where = — is the separation vector of the particles, with distance 
Tij = \Tij\, and unit vector = r^/r^. 

The dissipative part of the force acts to equalize velocities of nearby particles. 
It is also central force: 



pdiss 

ij 
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(5) 



The random part of the force represents a coupling to a heat bath: 



prand 



aw 



rand , 



Tij) dj ^ 



(6) 



where ^ = is a Gaussian random variable, independent for each ij pair of 
particles and timestep, with zero mean and At -1 variance. 



The coefficients 7 and a and the two weight functions cannot be chosen arbi- 
trarily: in order that the fluctuation-dissipation theorem holds [9], they must 
be related by 



w diss (r) 



.rand , 



W 
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(7) 



ct 2 = 2 7 A; B T . 



(8) 



We use the following weight function [10]: 



w diss (r) = w Tand (r) 



rand / 



\/l~r/r c , 
0, 



r < r n 



r > r n 



(9) 
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Given the interparticle forces, Newton's equations of motion are solved with 
a version of the velo city- Ver let algorithm [11]. As customary, our units were 
the cutoff length r c and the mass of a particle. 

In our plane Couette flow geometry (see Fig. 1) we have periodic boundary 
conditions in the x (streamwise) and y (spanwise) direction, and no-slip walls 
perpendicular to the z (gradient) axis. The walls are implemented as a soft 
repulsion potential in the normal direction: 



f 



wall 



-a wall z, z < 

0, < z < L (in bulk) , (10) 

a wall (z-L), z>L 



where L is the distance between the walls. To realize no-slip boundary condi- 
tions, at each timestep we update the particles' velocity component parallel 
to the walls as 

V|| <= V || +a(^)(v wall -V||), (11) 

where the factor a(z) is selected so that the velocity is unaffected in the bulk, 
the wall velocity is imposed upon particles well outside the walls, and there is 
a continuous crossover between these limits. Near the bottom wall it is 



a(z) = < 



0, z > (in bulk) 

—z/d , —d <z<0 (near wall) (12) 

1, z < —do (well outside wall) 



and similarly defined near the top wall. This approach is admittedly not ele- 
gant, as for example it has a non-trivial timestep dependence, but nevertheless 
achieves no-slip boundary conditions with minimal effort. 

With these boundary conditions the particles wander outside the nominal wall 
position to a limited distance depending on the temperature and pressure. For 
too sharp walls the density near the walls displays oscillations as a function 
of distance from the wall. These are probably the consequence of local crys- 
tallization near a sharp surface. This unphysical phenomenon is avoided by 
softening the wall; we used a wa11 = a. 

The z components of the stress tensor at the walls are readily available from 
the momentum transfer between the wall and the particles. In the bulk, how- 
ever, they have to be computed from the pairwise forces [12]: 

'<> = T7 { ~ Yl f U> r ^/3 - Yl - v) a (Vi - v)f, } . (13) 

! ' (id) 
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Fig. 2. (a) Viscosity as a function of the shear rate. For large shear rates the 
viscosity decreases slightly (shear thinning), (b) First normal stress difference and 
(c) first normal stress coefficient. The normal stress difference is quadratic for small 
shear rate 7, but the coefficient decreases significantly for 7 > 1. 

We validated our code for Newtonian fluids by calculating the stresses with 
different methods: at the walls, the global bulk value (where v is taken as the 
average velocity of particles at the same height z), and calculated locally [here 
Eq. (13) is integrated with a coarse graining kernel, which vanishes outside 
the neighborhood of a point]. All measurements were equal within error except 
the local stress, which differed by 10% — we attribute this to the difference in 
the contribution of fluctuations. We then calculated the viscosity from the 
shear stress, which compared well with an independent measurement from the 
relaxation time of the shear velocity profile's build-up at sudden start-up of 
the shear boundary conditions from rest. 



3 DPD for viscoelastic fluids 



To simulate viscoelastic hydrodynamics we connected DPD particles with 
springs: some fraction (in this paper all) of the particles are paired up to 
form dumbbells, so Eq. (3) receives an extra term f dumb . Linear (Hookean) 
springs give rise to nonphysical effects in elongational flow [2], therefore we 
used finite extensible non-linear elastic (FENE) springs: 

Mr 

U 1 - (r/r max ) 2 1 } 

where r is the separation between the endpoints of the dumbbell, H is the 
Hookean spring stiffness (at zero extension), and r max is the maximum exten- 
sion of the spring. 

As the first validation test, we measured the velocity profile in our shear cell, 
and obtained the expected linear profile. The shear viscosity, however, is now 
weakly dependent on the shear rate, as shown on Fig. 2a. 
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Fig. 3. Distribution of dumbbell configurations in the shear plane. Each dot cor- 
responds to the end-to-end vector of a dumbbell in the xz plane. For small 7 the 
orientation is isotropic, at moderate shear rates it is an elongated ellipse, and at 
large shear rates it is further deformed. The sides of the box represent the maximum 
elongation r max of the FENE springs. 

We also measured the first normal stress difference, see Fig. 2b. The figure 
shows that the model does achieve its goal, namely that it yields a normal 
stress difference, which for small 7 increases as 7 2 . At shear rate 7 > 1, the 
coefficient ^1 starts to decrease notably. 

The viscoelastic modes are represented by the FENE dumbbells, to which we 
have full access in the numerical simulations. Figure 3 shows the distribution 
of dumbbell orientation and extension. At low shear rates the dumbbells are 
isotropically oriented. At larger shear rates the dumbbells become elongated, 
with the typical orientation at a small angle with the stream direction. At the 
highest shear rate the distribution shows almost zero angle of the dumbbells 
with respect to the stream lines. 

The time evolution of a typical dumbbell orientation is plotted on Fig. 4. For 
small shear rates the trajectory resembles a random walk in the configuration 
space, while at larger shear rates the dumbbells tumble: if by fluctuation the 
two endpoints get to different height z, they are dragged with different ve- 
locities, resulting in a horizontal stretch. If at this point the endpoints move 
to the same height, the differential drag disappears, and the dumbbells can 
relax. This is augmented with the rotational component of the shear flow to 
yield a tumbling motion. 



4 Summary 

In conclusion, we presented preliminary results on extending DPD for vis- 
coelastic fluids by connecting DPD particles to form dumbbells. The model 
shows (shear rate dependent) normal stress difference, and a small amount 
of shear thinning. This is our first step to study numerically the instabilities 
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Fig. 4. Trajectories of the configuration of a dumbbell at different shear rates. A 
selected dumbbell is traced for 50 time units for the two smaller shear rates, and 20 
time units for the largest shear rate. 

and turbulence in viscoelastic fluids with a method complementary to direct 
numerical simulations of the constitutive equations. 
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